/** Author: David Powell
This file produces Figure 3 + Figure 6A + Figure A12 + Figure A13a	
**/


clear all
set more off
set mat 11000
set seed 8721

global dir "/jules/b/dpowell"
global DATA "${dir}/purdue/replication/DATA"
global OUTPUT "${dir}/purdue/replication/output"

use ${DATA}/finaldata


**convert to MEDs
gen oxycontin_rate=1.5*1000*oxycontin/(60*totpop) if year>1999 & year<2017
gen oxycodone_rate=1.5*1000*oxycodone/(60*totpop)
gen hydrocodone_rate=1000*hydrocodone/(60*totpop)


keep if year>1996

preserve
gcollapse *_rate [pw=totpop], by(year) fast
li
restore


preserve
gcollapse *_rate [pw=totpop], by(nontripl year) fast

set scheme s2color  

# delimit ;

twoway (connected oxycontin_rate year if nontripl==0 & year>1999 & year<2017, msize(vtiny) lpattern(dash) color(blue)  lwidth(thick) sort yaxis(1)) (connected oxycontin_rate year if nontripl==1 & year>1999 & year<2017, msize(vtiny)  lwidth(thick) color(red) sort yaxis(1))
, graphregion(color(white))   ylabel(,nogrid) xlabel(2000(2)2016, angle(45)) legend(label(1 "Triplicate States") label(2 "Non-Triplicate States") rows(1) order(2 1)) xtitle("Year") ytitle("Morphine Equivalent Doses Per Capita",  axis(1));
gr export ${OUTPUT}/fig3a.eps, replace;


twoway (connected oxycodone_rate year if nontripl==0 , msize(vtiny) color(blue) lpattern(dash) lwidth(thick) sort yaxis(1)) 
(connected oxycodone_rate year if nontripl==1, msize(vtiny) color(red)  lwidth(thick)  sort yaxis(1))(connected hydrocodone_rate year if nontripl==0 , msize(vtiny) lwidth(thick) sort yaxis(1) lpattern("___..") ) (connected hydrocodone_rate year if nontripl==1, msize(vtiny) lwidth(thick) sort yaxis(1) lpattern(longdash_dot) )
, graphregion(color(white))  xlabel(1997(2)2017, angle(45)) legend(label(1 "Oxycodone in Triplicate States") label(2 "Oxycodone in Non-Triplicate States") label(3 "Hydrocodone in Triplicate States") label(4 "Hydrocone in Non-Triplicate States") rows(4) order(2 1 4 3))  ylabel(,nogrid) xtitle("Year") ytitle("Morphine Equivalent Doses Per Capita",  axis(1));
gr export ${OUTPUT}/fig3c.eps, replace;


# delimit cr
restore



preserve
gen former_tripl=(stf==18 | stf==26)
replace nontripl=0 if former==1
gcollapse (median) oxycontin_rate , by(nontripl former year) fast

# delimit ;

# delimit ;
twoway (connected oxycontin_rate year if nontripl==1  & year>1999 & year<2017, lcolor(green) lpattern(dash_dot) lwidth(thick) msize(vtiny) sort yaxis(1))  
(connected oxycontin_rate year if former==1   & year>1999 & year<2017,  lwidth(thick)  lcolor(red)  msize(vtiny) sort yaxis(1))
(connected oxycontin_rate year if former==0 & nontripl==0  & year>1999 & year<2017,  lwidth(thick) lcolor(blue) lpattern(dash) msize(vtiny) sort yaxis(1))
, graphregion(color(white))  xlabel(2000(2)2016, angle(45))  legend(label(2 "Former Triplicate States") label(1 "Never Triplicate States") label(3 "Triplicate States")  rows(3)) xtitle("Year") ytitle("Median Per Capita OxyContin Supply" "Morphine Equivalent Doses" ,  axis(1)) ylabel(,nogrid);

gr export ${OUTPUT}/figA12.eps, replace;

# delimit cr
restore



***DIFFERENCES WITH CONFIDENCE INTERVALS

foreach nnn of numlist 1997/2017 {
	gen tt_`nnn'=nontripl*(year==`nnn')
}
tab year, gen(ttt) nof

local nnn=1
foreach var of varlist oxycodone_rate hydrocodone_rate {

	gen beta`nnn'=.
	gen se`nnn'=0	
	gen low`nnn'=.
	gen high`nnn'=.
	reg `var' tt_* ttt*  if year<2018 [aw=totpop], cluster(stf) 
	foreach nnn2 of numlist 1997/2017 {
		replace beta`nnn'=_b[tt_`nnn2'] if year==`nnn2'
		replace se`nnn'=_se[tt_`nnn2'] if year==`nnn2'
		boottest tt_`nnn2', boottype(wild) weighttype(webb) reps(9999) seed(23857389)
		matrix A=r(CI)
		replace low`nnn'=A[1,1] if year==`nnn2'
		replace high`nnn'=A[1,2] if year==`nnn2'		
	}
	local nnn=`nnn'+1
}

preserve
bys year: keep if _n==1
gen year2=year+.15
# delimit ;
twoway 
(rcap low2 high2 year2 if year>1996 & year<2018, color(green)  )  (connected beta2 year2 if year>1996 & year<2018, color(green) msize(vtiny) sort yaxis(1)) 
(rcap low1 high1 year if year>1996 & year<2018, color(red)  )  (connected beta1 year if year>1996 & year<2018, lpattern(dash) color(red) msize(vtiny) sort yaxis(1))
, graphregion(color(white)) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1997(2)2017, angle(45))   xtitle("Year") legend(off) ylabel(-2(2)8,nogrid) ytitle("Difference in Per Capita" "Morphine Equivalent Doses" ,  axis(1)) 
text(5 2014 "Oxycodone", place(e))
text(-1 2014 "Hydrocodone", place(e));
gr export ${OUTPUT}/fig3d.eps, replace;


# delimit cr
restore

gen beta`nnn'=.
gen se`nnn'=0	
gen low`nnn'=.
gen high`nnn'=.
reg oxycontin_rate tt_* ttt*  [aw=totpop], cluster(stf) 
foreach nnn2 of numlist 2000/2016 {
	replace beta`nnn'=_b[tt_`nnn2'] if year==`nnn2'
	replace se`nnn'=_se[tt_`nnn2'] if year==`nnn2'
	boottest tt_`nnn2', boottype(wild) weighttype(webb) reps(9999) seed(23857389)
	matrix A=r(CI)
	replace low`nnn'=A[1,1] if year==`nnn2'
	replace high`nnn'=A[1,2] if year==`nnn2'	
}
preserve
bys year: keep if _n==1
# delimit ;
twoway(rcap low3 high3 year if year>1999 & year<2017, color(red))  (connected beta3 year if year>1999 & year<2017, color(red) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) yline(0, lpattern(dash) lcolor(gs6)) xlabel(2000(2)2016, angle(45))  legend(off)  ylabel(,nogrid) xtitle("Year") ytitle("Difference in Per Capita OxyContin" "Morphine Equivalent Doses"  ,  axis(1));
gr export ${OUTPUT}/fig3b.eps, replace;


# delimit cr
restore


****FIGURE A13****
local nnn=`nnn'+1
gen beta`nnn'=.
gen se`nnn'=0	
gen low`nnn'=.
gen high`nnn'=.
reg oxycontin_rate tt_* ttt*  [aw=totpop] if nontripl==0 | stf==21 | stf==26 | stf==18 | stf==53 | stf==46, cluster(stf) 
foreach nnn2 of numlist 2000/2016 {
	replace beta`nnn'=_b[tt_`nnn2'] if year==`nnn2'
	replace se`nnn'=_se[tt_`nnn2'] if year==`nnn2'
	boottest tt_`nnn2', boottype(wild) weighttype(webb) reps(9999) seed(23857389)
	matrix A=r(CI)
	replace low`nnn'=A[1,1] if year==`nnn2'
	replace high`nnn'=A[1,2] if year==`nnn2'	
}
preserve
bys year: keep if _n==1
# delimit ;
twoway(rcap low4 high4 year if year>1999 & year<2017, color(red))  (connected beta4 year if year>1999 & year<2017, color(red) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) yline(0, lpattern(dash) lcolor(gs6)) xlabel(2000(2)2016, angle(45))  legend(off)  ylabel(,nogrid) xtitle("Year") ytitle("Difference in Per Capita OxyContin" "Morphine Equivalent Doses"  ,  axis(1)) ;
gr export ${OUTPUT}/figA13a.eps, replace;
# delimit cr
restore


****FIGURE 6***
drop tt_*
replace nontripl=0 if triplold==1
foreach nnn of numlist 2000/2016 {
	gen tt2_`nnn'=triplold*(year==`nnn')
	gen tt_`nnn'=nontripl*(year==`nnn')
}

gen beta=.
gen b2eta=.
reg oxycontin_rate tt_* tt2_* ttt*   [aw=totpop], cluster(stf) 
foreach nnn2 of numlist 2000/2016 {
	replace beta=_b[tt_`nnn2'] if year==`nnn2'
	replace b2eta=_b[tt2_`nnn2'] if year==`nnn2'
}


preserve
bys year: keep if _n==1
# delimit ;
twoway (connected beta year if year>1999 & year<2017, color(red)  lwidth(thick) msize(vtiny) sort yaxis(1)) 
 (connected b2eta year if year>1999 & year<2017, color(green)  lwidth(thick) lpattern(shortdash_dot) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) yline(0, lpattern(dash) lcolor(gs6)) xlabel(2000(2)2016, angle(45))  legend(label(2 "Former Triplicates") label(1 "Never Triplicates")  rows(1) order(2 1))  ylabel(0(.2)1,nogrid) xtitle("Year") ytitle("Difference in Per Capita OxyContin" "Morphine Equivalent Doses"  ,  axis(1));
gr export ${OUTPUT}/fig6a.eps, replace;


# delimit cr
restore

